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Abstract 

In this paper, we report the implementation of first-principles calculations of topological invariants 
Z2 within the full-potential linearized augmented plane-wave (FP-LAPW) formalism. In systems 
with both time-reversal and spatial inversion symmetry (centrosymmetric) , one can use the parity 
analysis of Bloch functions at time-reversal invariant momenta to determine the Z2 invariants. In 
systems without spatial inversion symmetry (noncentrosymmetric) , however, a more complex and 
systematic method in terms of the Berry gauge potential and the Berry curvature is required to 
identify the band topology. We show in detail how both methods are implemented in FP-LAPW 
formalism and applied to several classes of materials including centrosymmetric compounds Bi2Se3 
and Sb2Se3 and noncentrosymmetric compounds LuPtBi, AUTIS2 and CdSnAs2. Our work provides 
an accurate and effective implementation of first-principles calculations to speed up the search of 
new topological insulators. 

PACS numbers: 71.15.-m, 71.20.-b, 71.70.-d, 73.20.At 



I. INTRODUCTION 

Recently, topological insulators (TIs) have attracted great attention in the fields of con- 
densed matter physics and materials science. Based on the noninteracting band theory, TIs 
have gapped bulk gap and time-reversal symmetry protected metallic helical surface (edge) 
states where spin and momentum are locked together .l^ These novel physical properties hold 
great promise in applications of spintronics and quantum computing^I and have stimulated 
both experimental and theoretical studies. Indeed, the field of TIs is expanding so rapidly 
and there have been several excellent review articles on it.^^^ Although many TIs including 
alloy^, binary compounds^H^, ternary compounds^^"^, and quaternary compounds^ have 
already been theoretically predicted and experimentally realized, real materials that can be 
used in practical engineering are still needed. Therefore searching for new TIs with a vari- 
ety of excellent physical properties has become a central task in this filed. To achieve this 
goal, one has to develop an accurate and effective method to distinguish TIs from normal 
insulators. 

There are several general methods to determine the band topology of an insulator: 

(i) Based on the idea of bulk-edge correspondence of TIs,'^'^^ one can calculate surface 
(edge) states for a given insulator and count the number of gapless modes across the Fermi 
level. An odd number of gapless modes implies a TI while an even number indicates a 
normal insulator. This is a straightforward but not efficient way because the surface state 
dispersion may depend on every detail of the surface, for example, grown directions, termi- 
nated chemical elements and surface reconstructions. In some materials, the topologically 
nontrivial and trivial surface states can coexist, which further complicates the identifica- 
tion of the bulk topological order. To make sure that the gapless modes are topologically 
protected, one has to vary surface crystal structures and see if gapless modes can survive. 
Furthermore, a huge amount of computational resources is required in first-principles surface 
calculations. 

(ii) It is possible to use adiabatic continuity and so-called band inversion mechanism 
to identify TIs.^^™^ ^^ * ^^ ' The adiabatic continuity can be realized by artificially changing 
some external parameter such as the spin-orbit coupling (SOC) strength or lattice constant. 
Suppose the unknown state is near some known topological trivial or nontrivial states in a 
parameter space. If one tunes the parameter and the band gap stays open until it reaches 



the known state, then by the principle of adiabatic continuity, these two states share the 
same topological classification. Otherwise, the unknown state and known state may have 
different topological classifications if the band gap closes. Obviously, many intermediate 
calculations are required, making it a very tedious work. Band inversion at high-symmetry 
points within the Brillouin zone (BZ), as an empirical rule, can also be used to reveal the 
band topology. Although this empirical rule is adapted in some materials, such as half- 
Heusler^^"^, chalcogenide^^l^, and chalcopyrite^ compounds, it is not an universal way for 
arbitrary systems. 

(iii) The most general and direct approach is to calculate Z2 topological invariants from 
the knowledge of Bloch band theory.'^^"^ For materials with both time-reversal and spatial 
inversion symmetry (centrosymmetric systems), the simple parity criterion developed by Fu 
and Kane^^ have been applied in a number of -^orks.'I^IlSEsEelZS Qn the other hand, if the 
spatial inversion symmetry is absent (noncentrosymmetric systems), one must resort to a 
more complex method to evaluate Z2 invariants.'^^ Within a tight-binding framework, Fukui 
and HatsugaP^ have developed an effective algorithm to compute Z2 invariants in terms of 
the Berry gauge potential and the Berry curvature^ associated with the Bloch functions 
(BFs). This method has already been implemented in our first-principles codes and suc- 
cessfully predicted three-dimensional (3D) TIs in ternary half-Heusler-^- and chalcopyrite^Sl 
compounds and two-dimensional (2D) quantum spin hall effect (QSHE) in Silicene. Re- 
cently there appears another method which is in the same spirit of Ref . E3 but employs the 
charge center of Wannier functions.'^^'^ 

In this work, we illustrate the detailed implementation of first-principles calculations of 
topological invariants Z^ in both centrosymmetric and noncentrosymmetric systems within 
the full-potential linearized augmented plane-wave (FP-LAPW) formalism. Although the 
latter method for noncentrosymmetric systems can be applied to centrosymmetric systems, 
the parity criterion for centrosymmetric systems is a simpler and quicker way to determine 
the band topology. For this reason, we here introduce both of these methods. It should be 
emphasized that our methods are standard post-process after ground state wavefunctions 
are obtained in self-consistent calculation, so the calculation of Z2 invariants becomes a 
routine task just like band structures and density of states. Additionally, we have already 
paralleled our first-principles codes to speed up the calculation. Our implementation of the 
calculation of Z2 invariants is expected to be an efficient way for searching new TIs. 



The paper is organized as follows. In Sec. II, we review the fundamental expression 
of BFs within FP-LAPW formalism and the construction of overlap matrix, and then give 
the detailed formalism for implementation of parity analysis in centrosymmetric systems 
and lattice calculation of Z2 invariants in noncentrosymmetric systems. In Sec. Ill, we 
take centrosymmetric compounds Bi2Se3 and Sb2Se3 and noncentrosymmetric compounds 
LuPtBi, AUTIS2 and CdSnAs2 to illustrate the efficiency of our methods. Finally, we give a 
brief summary of our work in Sec. IV. In App. A, we provide details on the overlap matrix 
and its derivatives. 

II. METHODS 

In this section, we start by reviewing the formalism of BFs within FP-LAPW formalism 
and the construction of overlap matrix,'32iSl] then illustrate the calculation of Z2 invariants 
in both centrosymmetric and noncentrosymmetric systems. The key is to calculate the 
eigenvalues of parity operator according to parity criteriorP' (in the former case) or the 
overlap matrices related to time-reversal operator^ (in the latter case). 

A. Bloch functions and overlap matrix 

In the case of SOC, we consider BFs with two components. 
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where '|" and J, refer to the up and down component of spin. The periodic part of BFs is 
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where T is the transpose operator. The electrons in 
a solid environment have two different behaviors: those that are far from the nuclei and 
"free"-like can be described by plane waves, and those that are close to nuclei and unaffected 
by other nuclei can be described by atomic like functions. Within FP-LAPW formalism, the 
space is divided into two regions: a sphere with radius R^ around each atom, often called 
the muffin-tin region and the remaining space is interstitial region. ESlHll As a result, the BFs 
of electrons are always divided into two parts. Plane waves are used to construct the BFs 
in interstitial region 
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where Q is unit cell volume, z'^j^ ■ is the expansion coefficient, a and n stand for spin and band 
index, k for fc-points wave vector, Kj for the j-th reciprocal-lattice vector, j for the loop 
index of every expansion term and up to a largest value by the condition |A; + Kj\ < K^ax, 
and Kmax for the cutoff vector. Within the muffin-tin region (suppose the a-th atom sphere 
with radius Ra), the BFs can be written as 
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where r"" = r — t" and r" is the position of atom a; Im is the angular momentum index; 
Yim is spherical harmonics. In above formulas, Mp" = uf (r", -E"|) and -up" = iif (r°, -E'";^) 
are the radial solutions of scalar-relativistic Schri'/^cedinger equation of atom a and their 



u 



i,2 



energy derivatives, both evaluated at energy E^^. The local orbit radial functions 
uf (^r"',E^2) ^^6 added to the u^'" and u'^'^ for semi-core states (when / = Iq) and aimed 
to increase the variational freedom of standard basis functions. The last radial functions 
u'^y2 — ^1 ( '"") -^ri/2 ) ' ^^ ^^^ radial solution of full-relativistic Dirac equation, is also added 
to the uf and tt^'" but only for 5pi/2 or 6pi/2 orbits in heavy elements.'^ This extended 
full-relativistic local orbit can improve the accuracy of second-variational step when taking 
account of SOC. The AJ^ and 5^'^" are the coefficients of LAPW basis set, and B^is zero 



when APW basis set is used. AJ''^ , Bf'^ , C['^ , and -Df'^ are the coefficients of local orbit 
basis set. These coefficients can be determined by imposing various boundary conditions at 
the muffin-tin boundaries .'^SHlil 

Considering a lattice division within BZ, the overlap matrix between k point and its 
nearest-neighbor k -\- b has the form 



Overlap matrix Mmn is a very useful quantity in many Berry-phase related calculationSj^SH 
and the detailed formulas for its calculations are demonstrated in Appendix. 

B. Parity criterion in centrosymmetric system 

For systems with spatial inversion symmetry, Z2 invariants can be obtained by parity 
analysis developed by Fu and Kane^^. In 3D system there are eight time-reversal invariant 
momenta (TRIM) in BZ, rj=(„^„2n3) = ^i^iGi + n2G2 + n-^G-^), where Gj are primitive 
reciprocal-lattice vectors with rij = 0, or 1. The Z2 invariants are determined by the quan- 
tities 

N 



n^2™(r.). (6) 
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Here, C,2m i^i) = =tl is the parity eigenvalue of the 2m-th. occupied energy band at TRIMs 
Fj, i.e. {^2m,r, \P\ ^2m,ri); where P is parity operator. Because of the Kramers degeneracy 
at TRIMs, the 2m-th and (^m-i)-th occupied bands have the same eigenvalues, i.e., C,2m = 
$,2m-i- In 3D system, there are four independent invariants z/q; {1^11^21^3), given h^^ 

(-ir = ri<5. (7) 

{-'^Y" = [[ Si=(mn2n3), (8) 

where uo is independent of the choice of primitive reciprocal-lattice vectors Gj while z/i, z/2, 
and z/3 are not. A nonzero uq indicates that the system is a strong topological insulator 
(STI). When i^q = 0, the systems are further classified according to ui, 1^2 , and z/3. The 
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systems with z^i,2,or 3 7^ are called weak topological insulators (WTI), while 0; (000) is 
normal insulator (NI). 

To obtain Z2 invariants, the basic job is to calculate the matrix elements of parity operator 
(^nfe (^) |-P| ^nfc (^)) with even band index n at eight TRIMs Fj. The parity operator P is 
defined as {/;t}, where / is an inverse matrix making r — t- —r and t is a translational 
vector. Since parity operation will not change spin component of BFs, then, 

(vl>„fc (r) \P\ ^^^ (r)) = (^1, (r) \P\ ^l^ (r)) + (^i, (r) \P\ ^i^ (r)) . (9) 

In the following, we take ('ip^i^ (r) \P\ ilr^j^ (r) y as an example and suppress the spin index 
from here. Suppose that Tpnk (^) = Pi^nk (^) and inversion center at |, then ipnk (| ~ ■*") = 
i^nk (I + ^)- I^ can be rewritten as ipnk (^) = i^nk it — r), and finally we have Pipnk (^) = 
V'nfc (^ ~ ■*")• The matrix elements of parity operator are divided into two parts 

{^^k (r) \P\ V;„M (r)) = ii^nk (r) \P\ i^nk (r)), + J] (^^ M \P\ Ca, (r))^^^ . (10) 

a 

The contribution of interstitial region is 

{i^nk (r) \P\ ^„, (r)), = ^Y1 <k,^^-k, f e-^(^+^»)-^e'('=+^^-)-(*-'-) A (r) d\ 

= ^ E <fc.^"'^.e^^'^''^^*A (2fc + 1^. + i^,) . (11) 

Here, A(r) is a step function with zero value in muffin-tin sphere and unit value in 
interstitial region and A(i^) is its Fourier transformation. While inside the muffin- 
tin region, the radial coefficients in Eq. ^ can be rewritten as a product of two 
parts, one of which depends on atomic positions and the other does not. For example, 
A^ (n, k) = ^ • r]n,im {k + Kj, Ra) e*(*'+-'^j)''^", where r" is the position of a-th atom and 
Vn,im (k + Kj,Ra) is independent of r". Therefore, 

PAf^ (n, k) = J2vn,im (k + K,, R^) e^C^+^Mt-r-) _ (12) 

j 
If atom a is operated by parity operator, it must overlap another equivalent atom /3 by 

translated integer numbers of primitive real-lattice, i.e. t — r" = Rh + r^ with R^ = 



/iiHi + /i2a2 + h^a.^, where hj is integer number and a^ is primitive real-lattice vector. Then 
above equation can be rewritten as 

P^L in, k) = J2vnM (k + Kj, Rp) ^i^^J^.)-' ^^-^n 

3 

= AL (n, fc) e^'=■^^ (13) 

and there are similar operations for PB^^ {n, k), PCj^ (n, k), and PDf^ {n, k). The spherical 
harmonics operated by parity operator is, PYim [r") = (— 1) Yi^ {r°')- Then we have 

Pi^:k{r) = e'"-'"' Yl [^?- (^' ^^ <i + ^L in, k) u^ + C^ {n, k) ul, + D^ {n, k) u'^,/,] {-!)' Y,^ (f ^ 

Im 

(14) 
Therefore, we can easily obtain Pi'^ki''") ^J using radial coefficients of atom /3, A^^ {n, k), 
^hnin,k), Ci^{n,k), and Di^(n,k), which have already been calculated. Finally, the 
calculation of {iJ^^ (r) \P\ ip'^/, (r)) is very similar to {ip^k (^) li'nk (^))) which can be found 
in Appendix. 

C. Lattice calculation of Z2 invariants in noncentrosymmetric system 

A nontrivial topological invariant Z2 can be interpreted as an obstruction to make the 
BFs smoothly defined over BZ under time-reversal constrains.'^SHSZlfjgj-g^ -^g present a lattice 
evaluation of the Z2 invariants in terms of the Berry gauge potential and Berry curvature 
associated with the BFs.^^' This method has been recently applied to our first-principles 
studies of ternary half-Heusler^ and chalcopyrite^ TIs and QSHE in Silicene thinfilm.'^Sl 

We first briefiy describe the formalism for a 2D system. It was shown by Fu and Kane^ 
that under the time-reversal constraint, the Z2 invariants can be written as 



dk-Aik)- / d'kT{k) 

dB+ Jb+ 



mod 2, (15) 



^ 27r 

where B'^ and dB^ represent half of BZ and its boundary (Fig. QJ). The central quantities 
are the Berry connection 

A=z5^(u„(fc)| Vfcn„(fc)) (16) 
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Figure 1: (Color online) Schematic drawing of lattice mesh in a two-dimensional Brillouin zone. 
Under the time-reversal constraint, only half of Brillouin zone B'^ is needed, which is denoted by 
shaded region. The thick lines indicate the boundary of iS"*", i.e., dB'^ ^ and the open arrows denote 
their directions. All fc-points are divided into three classes: B'^ , Bj , and B'^, which are represented 
by small (black) solid, small (black) open and large (blue) solid circles, respectively. 

and the Berry curvature 

^(fc) = VfeX A(fc) U (17) 

where |-u„ (fc)) is the periodic part of BFs and the sum is over occupied bands. TIs are 
characterized hj Z2 = 1 while normal insulators have Z2 = 0. 

In the following, we introduce the calculation of \un (k)) in half of BZ referred to as B~^ 
([— Gi/2, Gi/2] (g) [— G2/2,0]) according to the time-reversal constraint. As shown in Fig. 
QJ the fc-points on a 2D BZ with N x N division are divided into three classes: Bf, B~ , and 
B'^. Firstly, we obtain |u„ (fc)) in B^ except for the points on the right edge. The points on 
the right edge (fc' = fc + Gi) are the periodic images of those on the left edge (fc), and can 
be calculated by using the periodic gaug#^'^ 



|m„ (fc + Gi)) = e-^^-'- |m„ (fc)) . (18) 

Secondly, we consider the Bj points on the boundary dB~^, i.e., the left part of the bottom 



edge and the right part of top edge. These points —k G B^ are the Kramers doublets of 
k G B~^ points, so they can be calculated by the time-reversal constraint, 

\un i-k)) = e \un (k)) , for fc G Bj. (19) 

where = —iayK is the time-reversal operator with K the complex conjugation. Note 
that translational phase factors must be properly considered. For example, k' G Bj and 
k G Bf are two points which are centrosymmetric about the midpoint of the bottom edge, 
i.e., k' = — fc — G2, then we have 

\Un{k')) = \u„i-k-G2)) 

= e*«-^ \un i-k)) 

= e'^^-^Q \un (fc)) . (20) 

Finally, we calculate |m„ (fc)) on TRIMs, i.e., i3°, satisfied by eH{k)e-^ = H{k). The 
eigenvalues are . . . e2n-i {k) = e2n (k) < e2n+i (k) = £2n+2 {k) . . . because of the Kramers 
degeneracy. In this situation, the time-reversal constraint is given by 

\u2n i-k)) = e \u2n-1 (k)) , -fc and fc G -B° • (21) 

There are six TRIMs in half of BZ B+, -Gi/2 - G2/2, -Gi/2, -G2/2, 0, Gi/2 - G2/2, 
and Gi/2. For the former four points, the 2n-th eigenstates can be obtained from (2?2— l)-th 
eigenstates by using above constraint. Here, one should also consider the translational phase 
factor, for example, 



\u2n{-Gj2-G2/2)) = e^(^^+^^)-^|M2„(Gi/2 + G2/2)) 

= e^(^^+^^)-^e \u2n-1 (-G1/2 - G2/2)) . (22) 

The other two points, Gi/2 — G2/2, and Gi/2, can be obtained by their periodic 
image points, i.e. |m„ (Gi/2 - G2/2)) = e"*^!'^ |m„ (-Gi/2 - G2/2)), |m„ (Gi/2)) = 
e-^^i-Mn„(-Gi/2)). 



After applying the time-reversal constrain Eq. (19) and Eq. (21) and periodic gauge Eq 



(18), we have obtained a new set of basis functions \un (k)). Next, we introduce the link 



variable that is central to many Berry-phase related calculations ,^2112 given by 
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U^ (k,) = N-' (kj) det {um (k,) I Un {kj + fi)) , (23) 

where A^„ -"^ (kj) = |det {um (kj) \ Un {kj + i-i))\ is the norniahzing factor and ^ is the unit vec- 
tor on the fc-mesh. In practice, {um (kj) \ Un {kj + fj,)) is the overlap matrix {um,k\un.k+ii) or 
its derivatives with the time-reversal operator B including {um,k\Qun,k+iJ.) , {Qum,k\un,k+iJ.) , 
and {Qum,k\Qun,k+ii) ■ The calculation of {um{kj) \ Un{kj + /j,)) is demonstrated in Ap- 
pendix. 

The finite element expressions for Berry connection .4. and Berry curvature J-" are 

A^{k,) = lm\ogU^{k,), (24) 

and 

J^ikj) = ImlogU^ (fc,) U^ {k, + ti) U^' (fc, + u) U,' (kj) , (25) 

where the return value of the complex logarithm function is confined to its principal branch 



(— 7r,7r]. We can then insert these expressions into Eq. (15) to calculate the Z2 invariants. 



To visualize the above procedure, an integer field n{kj) can be defined for each torus: 

nikj) = ^ {[A^A^ (fc,) - A^A. (fc,)] -J'ikj)} , (26) 

where A^ is the forward difference operator. The Z2 invariants are given by the sum of the 
n-field in half of the BZ, i.e., Z2 = Ylk gb+ ^i^j) niod 2. The sum of n-field configuration 
over the entire BZ gives a vanished Chern number for time-reversal invariant systems. It 
must be emphasized that the n-field summed over half of BZ is gauge-invariant module 2 
even though itself depends on a specific gauge choice. 

In 3D system, there are six possible 2D tori. These 2D tori are defined as follows: for 
example, the torus T{Xq) is spanned by G2 and G3 with the first component fixed at 0, and 
T{Xi) is obtained by fixing the first component at —Gi/2. The other four tori T{Yo), T{Yi), 
T{Zq), and T{Zi) are defined similarly. For each torus, one can calculate the corresponding 
Z2 invariants, xq, Xi, yo, yi, zq, and Zi, by using the steps outlined above for 2D BZ. Out 
of the six possible Z2 invariant only four of them are independent due to the constraint 
xq + xi = yo + yi = zq + z\ (mod 2). Following Refs. [3014321 we denote four independent 
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Figure 2: Schematic drawing of four independent tori in a three-dimensional Brillouin zone. The 
four independent tori T{Zq), T(Zi), T{Xq) and T(Yq) are located at k^, = 0, k^ = —G3/2, ki = 0, 
and /c2 = 0, respectively. 

Z2 invariants by uq; {viV2V?)i with z/q = {zq + Zi) mod 2, Vi = Xi, V2 = Vi and 2/3 = Zi. The 
corresponding four independent tori T{Zo), T{Zi), T(Xo) and T(Yo) are shown in Fig. 12} 



III. RESULTS 

In this section, we apply our methods to both centrosymmetric and noncentrosymmetric 
systems. In the case of centrosymmetric compounds Bi2Se3 and Sb2Se3, our parity analysis 
shows that Bi2Se3 is a STI while Sb2Se3 is a NI. The lattice calculation of Z2 invariants has 
also been used as a double check and the results are consistent with the parity analysis. We 
then turn to noncentrosymmetric compounds LuPtBi, AUTIS2 and CdSnAs2. By turning 
lattice constant, we studied three different topological phases of LuPtBi, i.e., STI, topological 
metal (TM), and NI. Furthermore, the Z2 invariants show that chalcopyrite compounds 
AUTIS2 and CdSnAs2 are STI and NI, respectively, in their native states without any strain. 

The calculations of band structures and Z2 invariants in this work were performed us- 
ing FP-LAPW method,'^^'^ implemented in the package wien2k.'^ We used two types 
of exchange-correlation potentials. The generalized gradient approximation (GGA)H^ was 
used for Bi2Se3 and Sb2Se3, while the modified Becke- Johnson exchange potential together 
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with local-density approximation for the correlation potential (MBJLDA)HS! was used for 
LuPtBi, AUTIS2, and CdSnAs2 because the resulting band topology is sensitive to the 
choice of exchange-correlation potentials in these systems.'^ The converged ground state 
was obtained using KmaxRMT = 9.0 for each system, where Kmax is the maximum size of 
reciprocal-lattice vector and Rmt represents the smallest muffin-tin radius. The fc-points 
sampling in BZ was also carefully checked such that self-consistent field calculations were 
well converged. Spin-orbit coupling was included by a second- variational procedure,'^ where 
states up to 9 Ry above Fermi level were included in the basis expansion, and the relativistic 
Pi/2 corrections^ were also considered for 5p and 6p orbit in order to improve the accuracy 
for systems including heavy elements. 

For a given system, the time taken by calculating of Z2 invariants depends on numbers of 
lattice divisions on four independent tori in 3D BZ and numbers of occupied bands considered 
below the Fermi level. For most of systems, a 10 x 10 lattice division on each torus is enough 
for obtaining a converged result just as mentioned in Ref. [Ml However, one must be very 
careful with the cases of small local band gaps, for example the system shown in Fig. [6tc), 
50 X 50 lattice division is need to reach the convergence. The included number of occupied 
bands should always been explicitly separated with other low-lying bands with an obvious 
global energy gap. The principle is that these low-lying bands are usually closed shell with 
much lower energy and should have trivial band topology. In the following, we chose 18, 18, 
30, 40 and 20 occupied bands for Bi2Se3, Sb2Se3, LuPtBi, AUTIS2 and CdSnAs2, respectively. 

A. Centrosymmetric systems 

To demonstrate the quality of our methods, we first test the centrosymmetric systems 
Bi2Se3 and Sb2Se3. Recently, Bi2Se3 family of compounds have been both theoretically and 
experimentally observed to be TIs with an exception of Sb2Se3.'^^'^ Tetradymite semicon- 
ductor Bi2Se3 family has a rhombohedral crystal structure with space group RSm (No. 166) 
and three nonequivalent atoms in a primitive cell. The calculated band structures of Bi2Se3 
and Sb2Se3 are presented in Fig. [3] with the lattice constants taken from previous studies.l^ 
The 18 occupied bands (— 6 ~ eV) are isolated from other low-lying bands and fully de- 
termine the topological nature of the systems, so we consider them as a bands group in the 
following calculation of Z2 invariants. 
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Bi2Se3 


-1 


+1 


+ 1 


+1 


+1 


+1 


+1 


+1 


Sb2Se3 


+ 1 


+1 


+ 1 


+1 


+1 


+1 


+1 


+1 



Table I: Parities 6i at eight TRIMs for Bi2Se3 and Sb2Se3. The relative coordinates in primitive 
reciprocal-lattice of eight TRIMs are (0, 0, 0), (0, 0, 0.5), (0, 0.5, 0), (0, 0.5, 0.5), (0.5, 0, 0), (0.5, 0, 
0.5), (0.5, 0.5, 0), (0.5, 0.5, 0.5). The Z2 invariants are 1; (000) for Bi2Se3 and 0; (000) for Sb2Se3, 
which indicate a STI and a NI respectively. 

Si §2 h 5i (55 ^6 ^7 <^8 1^0; (1^11^21^3) 

1;(000) 
0;(000) 

Because the existence of spatial inversion symmetry, the parity criteriorP^ is apphcable 
here. As a first step, we choose eight TRIMs in 3D BZ with relative coordinates (0, 0, 0), 
(0, 0, 0.5), (0, 0.5, 0), (0, 0.5, 0.5), (0.5, 0, 0), (0.5, 0, 0.5), (0.5, 0.5, 0), (0.5, 0.5, 0.5) 
in a primitive reciprocal-lattice. Then, we calculate the parity eigenvalues of 9 occupied 
bands with even band index (sorted by energy) out of 18 occupied bands at every TRIM. 
The parity of each TRIM, (5j=i^2,...8 iii Eq. (|6]), are obtained by multiplying over the parity 
eigenvalues of these 9 bands. The Z2 invariant z/q is obtained by multiplying over the parities 
of all TRIMs according to Eq. ([7^, while 2/^=1,2,3 by multiplying over the parities of TRIMs 
resided in the same plane according to Eq. (IS]). The 5i and Z2 invariants are listed in Table 
[ij The Z2 invariants are 1; (000) for Bi2Se3 and 0; (000) for Sb2Se3, indicating a STI and a 
NI respectively. One can see that the main difference lies at F point, i.e., 5i is —1 for Bi2Se3 
and +1 for Sb2Se3, while the other TRIMs share the same parities. We also give the parity 
eigenvalues of these 9 bands at F point, as listed in Table [IT| 

We have also used the lattice calculation of Z2 invariants as a double check. Figure |4] 
shows the n-field configuration for Bi2Se3. The Z2 invariants on each torus are zq = 1, Zi 
= 0, Xo = 1, and y^ = Ihy the sum of the ra-field in half of 2D BZ and then moduling 2. 
Total Z2 invariants 1;(000) indicate that Bi2Se3 is a STI. On the other hand, Figure [5] shows 
the ra-field configuration for Sb2Se3 with Zq = 0, Zi = 0, Xq = 0, and yo = on each torus. 
Total Z2 invariants 0;(000) indicate that Sb2Se3 is a NI. As expected, our lattice calculation 
of Z2 invariants are the same as parity analysis, and all of these two methods are consistent 
with the previous work.'^ 
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Figure 3: Band structures of strong topological insulator Bi2Se3 with Z2 invariants 1; (000) and 
normal insulator Sb2Se3 with Z2 invariants 0; (000). The eighteen occupied bands (every two of 
them are twofold degenerate) from —6 to eV are used to calculate Z2 invariants. The high- 
symmetry points in Brillouin zone are the same as Ref. 



Table II: Parity eigenvalues of Bi2Se3 and Sb2Se3 at V point for 9 occupied bands. The corresponding 
band energy increases from left to right. The parity of F point, (5i, is —1 for Bi2Se3 and +1 for 
Sb2Se3 respectively. 

5i 



Bi2Se3 
Sb2Se3 



-1 
-1 



-1 
-1 



-1 



-1 
-1 



-1 
-1 



B. Noncentrosymmetric systems 

Having established the effectiveness of our methods in centrosymmetric systems, we now 
turn to noncentrosymmetric systems by taking LuPtBi as the first example. It has been 
predicted that LuPtBi, as a member of ternary half-Heusler family, can realize a topological 
nontrivial state under uniaxial strain.'ISH^^E^Ml ^j^e crystal structure of LuPtBi is described 
by space group FAZm (No. 216) with three nonequivalent atoms in a primitive cell. The 
calculations were performed using the experimental lattice constant of 6.574 A."^ As shown 
in Fig. |6Va), LuPtBi is a semi-metal with small electron and hole pockets around Fermi level 
at F point. The band gap around F point can be obtained by applying an uniaxial strain, 
then 30 occupied bands (from —8 to about eV) were used to calculate Z^ invariants. 
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Figure 4: The n- field configuration for Bi2Se3 computed under tlie time-reversal constraints. The 
four tori are T{Zq), T{Zi), T{Xo) and T{Yq) with the shaded area indicating half of the 2D BZ. 
The white and black circles denote n = 1 and —1, respectively, while the blank denotes 0. The Z2 
invariants for each individual torus is obtained by summing the n-field over half of the torus and 
then moduling 2. These read zq = 1, ^i = 0, xq = 1, and yo = 1- The Z2 invariants of the system 
are 1;(000). 



As mentioned in our previous works,'^^ -^' topological phases of half-Heusler family are very 
sensitive to the change of lattice constants. Generally speaking, hydrostatic expansion leads 
to topological nontrivial phases while hydrostatic compression leads to topological trivial 
phases. Additionally, one must apply an uniaxial strain based on hydrostatic strain, i.e., a 
non-hydrostatic strain, to realize true topological insulating state because the states around 
Fermi level at F point are fourfold degenerate and protected by cubic symmetry. Therefore 
it is necessary to fully understand how the strain (hydrostatic and non-hydrostatic) acts on 
the topological phase in half-Heusler family. 

By turning lattice constants a{= h) and c, we found three different topological phases of 
LuPtBi including STI, TM, and NI, as shown in Fig. [6tb),[6tc), and|6](d), respectively. The 
non-hydrostatic strains can separate the fourfold degenerate states of valence and conduction 
bands around F point. In the case of Fig. [6tb), the global band gap together with Z2 
invariants 1; (000) indicate that this is a STI. While in the case of Fig. p[c), it is essentially 
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Figure 5: The n- field configuration of Sb2Se3. Tlie labels are the same as Fig. |4] The Z2 invariants 
for each individual torus read zq — 0, zi — 0, xq = 0, and yo — 0. The Z2 invariants of the system 
are 0;(000). 



a metallic state but has local band gap everywhere in the BZ. The Z2 invariants 1; (000) 
show a nontrivial state which is usually called TM. On the other hand, hydrostatic strain 
(large enough compression) can also create a band gap, just like Fig. [6td), but this is a NI 
because the Z2 invariants are 0; (000). 

Ternary chalcopyrite compounds of composition I-III-VI2 or II-IV-V2 are another impor- 
tant class of noncentrosymmetric TIs. In our previous work,!^ we have shown that a large 
number of ternary chalcopyrite compounds can realize the topological insulating phase in 
their native states. Here we take AUTIS2 and CdSnAs2 as noncentrosymmetric examples to 
show our methods for Z2 invariants calculation. The crystal structure of chalcopyrite is de- 
scribed by the space group IA2d (No. 122) with three nonequivalent atoms in a primitive cell, 
which can be regarded as a superlattice of the zinc-blende structure with small structural 
distortions. The crystal structure parameters of AUTIS2 rj = 1.016 and 6u = —0.018 are 
obtained by first-principles total energy minimization, and the experimental data rj = 0.980 
and 5u = 0.261'^ are used for CdSnAs2, where i] = c/2a is the tetragonal distortion ratio 
and Su is the internal displacement of anion.!^ AUTIS2 and CdSnAs2 are all semiconductors 
with band gap of 0.14 eV and 0.13 eV, as shown in Fig. [Wa) andjWb) respectively. Totally 
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(a) static lattice constant 



(b) a -4%a„ ; c -6%c, 



(I ' n 




Figure 6: Band structures of LuPtBi with the static lattice constant (a) [oq = 6o = cq = 6.574A], 
non- hydrostatic strains (b) [oq — 4%ao, cq — 6%Co] and (c) [ag + 6%ao, cq — 2%co], and hydrostatic 
strain (d) [ao — 8%ao, cq — 8%co]. The topological phases in (a), (b), and (c) are topological 
insulator, topological metal, and normal insulator, respectively. The 30 occupied bands (from —8 
to about eV) are used to calculate Z2 invariants. 

40 and 20 occupied bands (—6 ~ eV) are used to calculate Z2 invariants for AUTIS2 and 
CdSnAs2, respectively. We find that AUTIS2 is a STI with the Z2 invariants 1;(000) while 
CdSnAs2 is a NI with the Z2 invariants 0;(000). 

IV. SUMMARY 

In summary, we have presented the implementation of first-principles calculations of 
topological invariants Z2 in both centrosymmetric and noncentrosymmetric systems within 
FP-LAPW formalism. Generally, one can use a lattice version of Z2 invariants to identify 
the band topology, though in centrosymmetric systems, a simple parity criterion is possible. 
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Figure 7: Band structures of strong topological insulator AUTIS2 with Z2 invariants 1; (000) and 
normal insulator CdSnAs2 with Z2 invariants 0; (000) . The occupied bands which range from 
—6 ~ eV are included for calculating Z^ invariants, i.e. 40 bands for AUTIS2 and 20 bands for 
CdSnAs2 respectively. The high-symmetry points in Brillouin zone are the same as Ref. [5TJ 

The n-field configuration depends on a specific gauge, but the resulting Z2 invariants are 
gauge- invariant. Our method has two merits: (i) the algorithm implemented in our FP- 
LAPW framework is not expensive and the first-principles code can be easily paralleled; (ii) 
it is designed as a standard post-process of first-principles calculations, so the identification 
of topological nature for a given material becomes a routine task. Therefore, our method is 
able to identify TIs in relatively short time and we anticipate it will speed up the discovery 
of new topological insulators in future. 

Acknowledgments 

This work is supported by NSF of China (Grants No. 10974231) and the MOST Project of 
China (Grant No. 2007CB925000 and 2011CBA00100). W.F. was supported by the LDRD 
Program of ORNL. D.X. acknowledges support by the Materials Sciences and Engineering 
Division, Office of Basic Energy Sciences, U.S. Department of Energy. We also acknowledge 
the computational resources supported by Texas Advanced Computing Center (TACC) and 
Supercomputing Center of Chinese Academy of Sciences (SCCAS). 



19 



Appendix A: Overlap matrix and its derivatives with the time-reversal operator 

In this appendix, we give the overlap matrix {um,k\un,k-irb) and its derivatives with the 
time-reversal operator 6, including {um,k\'S)Un,k+b) , {Qum,k\un,k+b) , and {<dUm,k\Qun,k+b) , 
where m and n stand for band indexes, and b stands for unit vector /x or i/ on fc-mesh (see 

Fig. ig. 

Firstly, we consider the overlap matrix {um,k\un,k+b) according to Eq. (IS]), 



{Um,k\Un,k+b) 



t I t 
'^m,k\'^n,k+b 



+ {uiJu^ 



m,k\"'n,k+b / ■ 



(Al) 



Here we only discuss (umkl^nk+b) because that ( w^ /, | ^„ fc_|_5 ) has the similar formulas. 
Like the BFs, the overlap matrix can also be divided into two parts: interstitial region and 
muffin-tin region. 



t I t 

'^m,k\'^n,k+b 



^m,k\'^n,k+b / + ^^ \^m,k\'^n,k+b 
a 

The contribution of interstitial region is 



MT° 



(A2) 



u, 



,/ t/ cell 

J 3 

= J2T.<'^A,k+b,r^iK,-K,,). 

J f 



(A3) 



Here, A(r) is a step function, it have zero value in muffin-tin sphere and unit value in 
interstitial region, and it's Fourier transform is 



a 

The contribution of a-th muffin-tin sphere is 



-iK-T-^T^RljliKRc 



VL KR„ 



t I t 

'^m,k\'^n,k+b 



MT'^ 



,ifc-(T"+r°) 



MT°' 
bT° 



^l^Ar) 



-J(fc+6)-(T"+r") /ta 



ifa 



^rfc(^) ^n 



Using the Rayleigh plane-wave expansion 



,,tQ 
i,k+b 



^!:k+bir)d'r 



(A4) 
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-ibr° 



A-nY. (-'f ^V'^" (^) ^'"™" (^")>' (^^") 



V'm" 



where r" = |r°|, 6 = |b|, and jin {br°') is the spherical bessel function. Then, 



(A5) 



t I t 



MT° 



X 



4vre-*-5^(-.)'"yir,„„(6 

l"m" 



Im I'm' 
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t« 
Zm 
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t« 
Zm 

fa 
Im 
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im 



D 



D 



D 



D 



{m, 


fc) 


{m, 


fc) 


(m, 


k) 


(m, 


k) 


(m, 


k) 


(m, 


k) 


(m, 


k) 


{m, 


k) 


(m, 


k) 


(m, 


k) 


(m, 


k) 


(?72 


k) 


{m 


k) 


{m 


k) 


(m 


k) 



Bll, {n,k + b 
C]L' in,k + b 
Dj'l' in,k + b 



4m' in,k + b 

Bll' in,k + b 

4m' in,k + b 

Dll' in,k + b 

^fw in,k + b 

Bll, {n,k + b 

Cll' in,k + b 

Dll, {n,k + b 

^ ^Jm' in,k + b 

^ BJZ' in,k + b 

^ Cll' in,k + b 

^ Dll' in,k + b 



n,k + b) 



«/ 1 ^'AJl",b 



f,a ■f,a ■ 
f.a f.a ■ 



^1,1 ^l',l/23l",b 
■ t.a t.a ■ 

• t,ct • f,a ■ 
Hi HAJl",b 






"i,l H, 1/231", b 



^i,2 «iMJ/",fe 
t,Q . t,a ■ 

t,o t," ■ 
«/,2 W/',2J«",b 



^1,2 H,l/23l",b 

t,a ti" ■ 
Hl/2H,l3l",b 

t,a . t,a ■ 
'^1,1/2^,131", b 

t,a t,a ■ 
^1,1/2^,231", b 



'^l,l/2''^l', 1/2^1", b 



Gmm 
ll'l" 



WV' 



(A6) 



Therefore, the matrix elements Iv} .l-ul 



''m.fc I "n.fc+b 



MT-^ 



are constructed by two parts: radial in- 
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tegrals and angular integrals. The radial integrals are 
























"«,1 ^l',l/23l",b 



«/,2 «;' ,lJ«",b 



«/ 2 «;' lJ«",b 
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2 t /^Q TTia 



r -u r 



;( 



i?" 



2„,t l^a 



r u [r 



f( 



i?" 



2 t 



4(^ 



r^-u[ (r" 



i?" 



r^-uj fr° 



4(^ 
•t /' 



r^-uT fr" 



i?" 



r^-uj ir'^ 



2.-,t /'^Q 



r M r 



n 



R" 



r'^uj [r" 



4(^ 



2.,t /'^« 



r M r 



H 



R°' 



.2,,t 



;(' 



r -u r 



2.,t /'^« 



r M r 



H 



i?" 



r^-uj (r" 



2„,t /^Q 



r u [r 



l( 



R°' 



2„,t f^Oi 



r u [r 



f( 



i?" 



2„,t /^Q 



r M r 



K 



Er,i 
Ki 

K,i 
K,i 

Ki 



^«,2 



rpa 



rpa 
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u],{t^,EI,)ji. 
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u 
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1{t-,EI,)j,, 
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{hr"") dr, 
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I' {br°') dr, 
{hr'^) dr, 
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Er,y,)4{r'^,El,)jiiiibr")dr, 



EZ,/,)ul{r-,El,/,)jiii{hr-)dr, (A7) 



and the angular integrals is the Gaunt coefficients 



G— - = j Y,l^ (r) Yvra' (r) Y.i^n (r) d^. 
Secondly, we consider the matrix element {um,k\Qun,k+b) , 



(AJ 
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(M™,fc|eMn,fc+6) = - \u\^Ml^+b) + \^L,kWn,k+b) ^ (A9) 

and take ( w„jfc|M„*;^,_|_f,\ as example, it can be divided into two parts: 



''m,k\'^n,k+bj ~ \^m,k\'^n,k+bj + /-^ V™''" ''^+^/mT« ' \-^^^) 

ex 

Within the interstitial region, 



I 4-* \ _ V^ V^ t* 4-* jL /" f- 

,,fcl'"?i,fc+6/ ~ / . / . ^mk,j^nk+b.j' o / ^ 
' ■' ■ ., ^' J cell 



T j' 

E E ^l*-..^iw'0 (^. + ^/) ' (All) 



while inside the muffin-tin region (a-th atom sphere), 
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Thirdly, we consider the matrix element (0'Um,fc|'"n,fc+6)) 



(0Wm,fc|w?i,fc+6) 



4-* I t \ I / t* I 4- 
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Within the interstitial region, 
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{n,k -\- b 
{n,k + b 
{n,k -\- b 
{n,k + b 
(n, fc + 6 
{n,k + b 
(n, fc + 6 
(n, fc + b 
{n,k + b 
{n,k + b 
(n, fc + b 
(n, fc + 6 



4-,o t,Q ■ 

Ui I Ui, iJl",2k+b 



'^1,1 '^l',lJl",2k+b 
""/ 1 ""f 2Jl",2k+b 






• 4-, a tiC* ■ 

""(,1 '^l',lJl",2k+b 

• i,a • t,Q ■ 

""«,! '^l',lJl",2k+b 

""/ 1 % 2Jl",2k+b 






"^1,2 '^l',lJl",2k+b 

l,a • t,Q ■ 
'"«,2 '^l',lJl",2k+b 

'"/,2 '^l'2jl",2k+b 



^1,2 Ui, 1/2^1", 2k+b 

l,a t,a ■ 
"^1,1/2^1', l3l",2k+b 

\.,ct • t,a ■ 
"^1,1/2^1', l3l",2k+b 

^l,l/2^l',2Jl",2k+b 

'^l,l/2'^l',l/23l",2k+b 



(A15) 



}(-irG- 



mmm 
I" 



(A16) 
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Finally, we consider the matrix element (O'Um,fc|0'Wn,fe+f)), 



and take \uj^ik\^nk+b/ ^^ example, it also can be divided into two parts 



^'"m,fcl^n,fc+6/ — \^m,k\'^n,k+b/ ^ + 2-^ \^m,k\'^n,k+b/ ^^^ • (^18) 

a 

Within the interstitial region. 



J j' 



while inside the muffin-tin region (a-th atom sphere), 
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■I* I 4-* 



„,.=*-*'"E'"'>'?-(s 



l"m" 



Im I'm' 






+ B] 



+ Bt 



+ D 



' Im 



+D 



Im 



m, k) 



^I'm' 



m 



m, k) [d|4, 

^I'm' 
^I'm' 
^I'm' 



m, k) 
m, k) 
m, k) 



m 



k) [A\r^, 
m, k) \Bfr^, 



m, k) 



ia 



^I'm' 



m, 



k) ^I'm' 



m, k) 



ia 



B^ 

^I'm' 



m, k) [C,t, 



m, k) 



^I'm' 



{n,k + b 
{n,k -\- b 
{n,k + b 
(n,k + b 
{n,k + b 
{n,k + b 
{n,k + b 
(n,k + b 
{n,k + b 
(n, fc + 6 
{n,k + b 
{n,k -\- b 
{n,k + b 
{n,k + b 
{n,k + b 



n 



^,k + b) 



i.a i,a ■ 



I.a • l,a ■ 

W,i ^',iJi",b 

I.a l,a ■ 
<1 Ul',2Jl",b 

4,, a 4,, a 
^1,1 ^l',l/23l":b 

' I.a I.a ■ 

"ti K,iJi",b 

• 1,0 I.a ■ 



■ 4-, a -1,0 

^1,1 ^l',l/23l",b 

I.a l,a ■ 
<2 Ul'\lJl",b 



^'2 ^',iJi",b 

I.a l,a ■ 
<2 Uj,\2Jl",b 



4, a -lyCe 
^1,2 ^l',l/23l",b 

\l/2^l',l3t",b 

I.a • 1,Q ■ 
^l,l/2^i',lJl",b 

4, a -lyCe ■ 
^1,1/2^1', 2Jl",b 

4, a 4i'* 
'^l,l/2'^l',l/23l",b 



f_1 \'"+™' f-^-m-m'm," 

(A20) 



* Electronic address: ygyao@iphy.ac.cn 

1 C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005). 

2 C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801(2005). 

3 J. E. Moore, Nature (London) 464, 194 (2010). 

^ X.-L. Qi and S.-C. Zhang, Physics Today 63, 33 (2010). 



5 X.-L. Qi and S.-C. Zhang arXiv: 1008.2026 



27 



^ M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010). 

■^ J. C. Y. Teo, Liang Fu, and C. L. Kane, Phys. Rev. B 78, 045426 (2008). 

^ D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature (London) 

452, 970 (2008). 
9 B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006). 
^^ M. Konig, S. Wiedmann, C. Briine, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and 

S.-C. Zhang, Science 318, 766 (2007). 

11 L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007). 

12 H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature Phys. 5, 438 (2009). 

13 Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, 
and M. Z. Hasan, Nature Phys. 5, 398 (2009). 

14 Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. 
Dai, Z. Fang, S. C. Zhang, L R. Fisher, Z. Hussain, and Z.-X. Shen, Science 325, 178(2009). 

15 D. Xiao, Y. Yao, W. Feng, J. Wen, W. Zhu, X.-Q. Chen, G. M. Stocks, and Z. Zhang, Phys. 
Rev. Lett. 105, 096404 (2010). 

16 S. Chadov, X. Qi, J. Kiibler, G. H. Fecher, C. Felser, and S. C. Zhang, Nature Mater. 9, 541 
(2010). 

1^ H. Lin, L. A. Wray, Y. Xia, S. Xu, S. Jia, R. J. Cava, A. Bansil, and M. Z. Hasan, Nature Mater. 

9, 546 (2010). 
1^ H. Lin, R. S. Markiewicz, L. A. Wray, L. Fu, M. Z. Hasan, and A. Bansil, Phys. Rev. Lett. 105, 

036404 (2010). 
1^ B. Yan, C.-X. Liu, H.-J. Zhang, C.-Y. Yam, X.-L. Qi, T. Frauenheini, and S.-C. Zhang, Europhys. 

Lett. 90, 37002 (2010). 

20 Y. L. Chen, Z. K. Liu, J. G. Analytis, J.-H. Chu, H. J. Zhang, B. H. Yan, S.-K. Mo, R. G. 
Moore, D. H. Lu, L R. Fisher, S. C. Zhang, Z. Hussain, and Z.-X. Shen, Phys. Rev. Lett. 105, 
266401 (2010). 

21 T. Sato, K. Segawa, H. Guo, K. Sugawara, S. Sounia, T. Takahashi, and Y. Ando, Phys. Rev. 
Lett. 105, 136802 (2010). 

22 K. Kuroda, M. Ye, A. Kimura, S. V. Eremeev, E. E. Krasovskii, E. V. Chulkov, Y. Ueda, K. 
Miyamoto, T. Okuda, K. Shimada, H. Namatame, and M. Taniguchi, Phys. Rev. Lett. 105, 
146801 (2010). 

28 



23 W. Feng, D. Xiao, J. Ding, and Y. Yao, Phys. Rev. Lett. 106, 016402 (2011). 

24 Y. Sun, X.-Q. Chen, S. Yunoki, D. Li, and Y. Li, Phys. Rev. Lett. 105, 216406 (2010). 

25 B. Yan, H.-J. Zhang, C.-X. Liu, X.-L. Qi, T. Frauenheim, and S.-C. Zhang, Phys. Rev. B 82, 
161108 (2010). 

26 J. Kim, J. Kim, and S.-H. Jhi, Phys. Rev. B 82, 201312 (2010). 

27 H. Jin, J.-H. Song, A. J. Freeman, M. G. Kanatzidis, Phys. Rev. B 83, 041202 (2011). 

28 H.-J. Zhang, S. Chadov, L. Miichler, B. Yan, X.-L. Qi, J. Kiibler, S.-C. Zhang, and C. Felser, 
Phys. Rev. Lett. 106, 156402 (2011). 

29 S. Chen, X. G. Gong, C.-G. Duan, Z.-Q. Zhu, J.-H. Chu, A. Walsh, Y.-G. Yao, J. Ma, and S.-H. 
Wei, Phys. Rev. B 83, 245202 (2011). 

30 L. Fu, C. L. Kane and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007). 

31 J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (2007). 

32 R. Roy, Phys. Rev. B 79, 195322 (2009). 

33 L. Fu and C. L. Kane, Phys. Rev. B 74, 195312 (2006). 

34 T. Fukui and Y. Hatsugai, J. Phys. Soc. Jpn. 76, 053702 (2007). 

35 D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010). 

36 C.-C. Liu, W. Feng, and Y. Yao, |arXi v:1104.12 90l Phys. Rev. Lett, (in press). 
3^^ A. A. Soluyanov and D. Vanderbilt, Phys. Rev. B 83, 235401 (2011). 

38 R. Yu, X.-L. Qi, A. Bernevig, Z. Fang, and X. Dai, |arXiv:1101.2011 

39 D. J. Singh, Planewaves, Pseudopotentials and the LAPW Method (Kluwer Academic, Boston, 
1994). 

40 S. Bliigel and G. Bihlmayer, John von Neumann Institute for Computing, NIC Series 31, 85 
(2006). 

41 P. Blaha, K. Schwarz, G. Madsen, D. Kvaniscka, and J. Luitz, Wien2k, An Augmented Plane 
Wave Plus Local Orbitals Program for Calculating Crystal Properties (Vienna University of 
Technology, Vienna, Austria, 2001). 

42 J. Kunes, P. Novak, R. Schmid, P. Blaha, and K. Schwarz, Phys. Rev. B 64, 153102 (2001). 

43 R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993). 

44 R. Resta, Rev. Mod. Phys. 66, 899 (1994). 

45 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 

46 F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009). 

29 



'^'^ W. Feng, D. Xiao, Y. Zhang, and Y. Yao, Phys. Rev. B 82, 235121 (2010). 

"^8 W. Al-Sawai, H. Lin, R. S. Markiewicz, L. A. Wray, Y. Xia, S.-Y. Xu, M. Z. Hasan, and A. 

Bansil, Phys. Rev. B 82, 125208 (2010). 
^9 M. G. Haase, T. Schmidt, C. G. Richter, H. Block and W. Jeitschko, J. Sohd State Chem. 168, 

18 (2002). 
^^ B. R. Pamplin, T. Kiyosawa, and K. Masumoto, Prog. Cryst. Growth Charact. 1, 331 (1979). 
^1 S. Limpijumnong, and W. R. L. Lambrecht, Phys. Rev. B 65, 165204 (2002). 



30 



